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ABSTRACT 


The  design  of  surfaces  of  vehicles  is  a  lengthy 
and  tedious  process.  Advances  in  interactive  computer 
graphics  and  mathematical  representation  of  curves  and 
surfaces  have  made  interactive  graphics  an  attractive 
medium  for  design.  The  Surface  and  Modeling  and 
Analysis  program  provides  a  set  of  computer  sub¬ 
routines  for  (1)  calculating  descriptive  information 
about  a  surface  and  (2)  manipulating  data  (by  providing 
the  interface  between  the  programmer  and  the  data 
structure) .  These  routines  are  based  upon  parametric 
cubics  and  the  Coons'  Patch  technique.  A  brief  de¬ 
scription  of  the  representation  is  included. 

The  SMA  subroutine  package  was  implemented  on 
the  CDC  6700  Computer. 


ADMINISTRATIVE  INFORMATION 

This  work  was  performed  within  the  Graphic  Systems  Development 
Group,  Computer  Sciences  Division,  under  Task  15324,  Subtask  Area 
SR0140301,  Job  Order  1-1802-001  the- Mathematical  Sciences  Subroutine 
Library  Project. 
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INTRODUCTION 


Traditionally,  surface  design  has  been  carried  out  on  the  drawing 
board.  With  recent  developments  in  the  area  of  interactive  computer 

graphics,  and  with  improvement  in  plotters  and  numerically  controlled 

*- 

machine  tools,  new  methods  of  computer-aided  geometrical  design  have 
emerged.  These  methods  are  based  on  mathematical  descriptions  of 
shape  from  which  drawings  or  instructions  for  a  numerically  controlled 
machine  tool  can  be  produced.  The  non-mathematical  representation  of 
a  geometric  shape — for  example,  by  a  grid  of  points — is  both  wasteful 
of  storage  space  in  the  computer  and  ambiguous,  since  a  complete  de¬ 
scription  of  the  object  is  not  available.  The  method  of  mathematical 
representation  selected  must  be  capable  of  describing  the  kinds  of 
shapes  which  are  encountered  in  design.  The  method  must  have  the 
ability  to  generate  rapidly  any  particular  view  of  a  surface,  and  it 
must  be  axis^independent.  Some  mathematical  representations  do  not 
permit  multiply-valued  curves  with  large  slopes  relative  to  the  co¬ 
ordinate  axis,  or  closed  curves.  Parametric  methods  of  curve  and 
surface  description  have  been  widely  used  to  overcome  these  restrictions. 
Curves  are  represented  by  vector-valued  functions  of  a  parameter,  t,  for 
example 

P(t)  =  [f  (t)  ,g(t)  ,h(t)] 

and  surfaces  are  represented  by  functions  of  two  parameters,  u  and  v, 
for  example 

P(u,v)  =  [f  (u,v)  ,g(.u,v)  ,h(u,v)] 
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In  general,  the  shape  of  a  real  object  is  non-analytic  so  that  local 
conditions  do  not  have  an  overall  effect  on  shape.  For  this  reason 
various  piece-wise  representations  of  curves  and  surfaces  have  been 
devised  by  which  the  extent  of  the  influence  of  a  local  change  in 
shape  can  be  controlled. 

A  mathematical  representation  can  be  obtained  in  either  of  two 
basic  ways:  by  fitting,  or  by  design.  In  the  former  method,  the 
problem  is  to  obtain  a  mathematical  representation  for  a  shape  that 
has  not  been  mathematically  defined  but  which  exists  either  as  a 
physical  model  from  which  data  points  can  be  measured  or  as  the 
result  of  some  procedure.  Since  the  shape  is  already  defined,  the 
fitting  process  can  be  largely  automatic.  In  the  design  method,  how¬ 
ever,  the  problem  is  either  to  create  a  shape  which  satisfies  design 
constraints,  or  to  modify  an  existing  mathematically  defined  shape. 

In  this  case,  human  intervention  is  often  essential. 

The  Surface  Modeling  and  Analysis  (SMA)  subroutine  package  uses 

12  3 

both  parametric  cubics  and  "Coons'  Patch"  Techniques  *  *  .  The  first 
version  of  SMA  (Modi,  described  here)  may  be  used  in  first-order  design 
of  a  surface  (mathematical  modeling) .  A  subsequent  version  (SMA-Mod2) 

^  Coons,  S.  A.,  "Surface  for  Computer-Aided  Design  of  Space  Forms," 

MIT  Project  MAC,  MAC-TR-41  (Jun  1967). 

2 

Forrest,  A.  R. ,  "Curves  and  Surfaces  for  Computer-Aided  Design," 
Cambridge  University,  CAD  Group,  Ph.D.  Thesis  (Jul  1968). 

3 

Armit,  A.  P.,  "A  Multipatch  Design  System  for  Coons'  Patch,"  IEE 
International  Conference  of  Computer-Aided  Design,  Conference 
Publication  51,  Southampton  (Apr  1969). 
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planned  for  development  later  in  1973  will  enable  the  user  to  design 
a  new  surface  under  certain  constraints  and  will  provide  a  simplified 
method  of  modification  of  an  existing  mathematical  model.  It  will  also 
provide  capabilities  for  analysis  of  the  model  (area,  moments,  etc.)  created 
by  SMA. 


SMA  SUBROUTINE  PACKAGE 

The  functions  performed  by  the  various  SMA  component  routines  are 
summarized  in  Table  1  on  the  opposite  page.  The  package  enables  (1) 
the  manipulation  of  data  (via  the  interface  it  provides  between  the 
programmer  and  the  data  structure) ,  and  (2)  calculation  of  descriptive 
information  about  the  surface  (interior  and  boundary  points,  tangents, 
and  twist  vectors).  At  present,  the  data  structure  resides  in  a  pro¬ 
grammer-defined  array;  later  versions  of  SMA  will  use  disk  files. 

A  typical  surface  and  an  enlargement  of  a  single  patch  are  il¬ 
lustrated  in  Figure  1,  page  6,  using  SMA  patch  notation. 

Calling  parameters  and  variables  used  extensively  in  the  SMA 
subroutines  are  explained  in  Table  2,  page  7. 


4 


TABLE  1  -  SUMMARY  OF  SMA  SUBROUTINE  CAPABILITIES 


INTERFACE  ROUTINES 

LGETV.LSETV 

Retrieves  (or  stores)  surface  patch  data 
from  (or  into)  the  array  PATCH  according 
to  given  patch  numbers. 

IGETV.ISETV 

Retrieves  (or  stores)  surface  patch  data 
from  (or  into)  the  array  PATCH  using 
given  curve  numbers  at  intersections. 

MAIN  COMPUTATIONAL  ROUTINES 

SPLINE 

Computes  endpoint  rates  (related  to 
tangents)  and  puts  them  into  the  data 
structure. 

GENSURF 

Generates  boundary  and  interior  points 
of  a  patch. 

GENERAL  PURPOSE  SUBROUTINES 

CUBIT 

Fits  a  parametric  cubic  spline  through 
a  series  of  points. 

ICONP 

Compare  3-dimensional  vectors  for 

ICONT 

proportionality  and  identicality. 

EPTC 

Checks  for  user-defined  tangents  and  sets 
up  an  endpoint  tangent  condition  array. 

MZCAM 

Compares  a  3-dimensional  vector  to 
(-0.,  -0.,  -0.) 

MULT 

Performs  matrix  multiplication. 

BNDRYPT 

Determines  points  on  the  boundary  of 
a  patch. 

FGBLEND 

Determines  the  Coons '  blending  function 
value  at  a  given  parameter  value. 
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PATCH  (N,i,j) 


NV 
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(Corner  3)^  /  V  =  1 


Patch  Notation: 


PATCH (N, IP, JP)  where  N  denotes  an  entire  patch 

IP,JP  describe  the  topo¬ 
logical  position  of  the 
patch 


Point  Notation: 


Tangent  Vectors : 


Twist  Vectors: 


where  i  is  the  U-curve  number 
j  is  the  V-curve  number 


Tuij  9u^P(u,v))  | u , v  =  0  or  1 

Tvij  •  l^(P(u-v))  |UfT  .  0  or  1 
32 

W,  .  =  -  -(p(u,v) )  n  n 

ij  3u3v  u,v  =  0  or  1 


Figure  1  -  SMA  Notation 
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TABLE  2  -  EXPLANATIONS  OF  SURFACE  PATCH  NOTATION 


Variable 

or 

Array 

Description 

MAX 

u 

The  number  of  curves  in  the  U-direction  to  be  contained 
dn  the  data  structure.  (See  Figure  1.) 

MAX 

V 

The  number  of  curves  in  the  V-direction  to  be  contained 
in  the  data  structure.  (See  Figure  1.) 

Ml 

MAX  -  1 
u 

M2 

MAX  -  1 

V 

PATCH 

The  name  of  the  array  containing  the  surface  patch  data 
structure.  The  array  is  dimensioned  PATCH(48,  Ml,  M2). 

(See  SMA  Data  Structure) . 

IP,  JP 

Subscripts  of  specific  patch  to  be  referenced:  (PATCH(n, 
IP,  JP)  where  n  denotes  a  vector  component. 

1  x  IP  i  Ml,  1  £  JP  £  M2.  When  referring  to  an  entire 
patch,  N  will  be  in  for  n. 

DATA  STRUCTURE  INTERFACE  ROUTINES 

Each  Coons'  patch  consists  of  16  vectors:  four  point-position 
vectors,  eight  tangent  vectors,  and  four  twist  vectors.  An  SMA  patch 
consists  of  48  sequential  locations  in  an  array.  The  patch  is  divided 
into  four  partitions,  one  for  each  corner  point.  Each  corner  point 
partition  is  subdivided  into  four  3-dimensional  vectors.  Within  the 
partition,  the  vectors  are  stored  in  the  following  order:  (1)  point; 
(2)  tangent  along  a  U-curve;  (3)  tangent  along  a  V-curve;  (4)  twist 
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vector.  The  components  of  the  vectors  can  be  stored  in  any  order, 
but  the  order  must  be  consistent  throughout  the  structure.  The 
entire  system  of  patches  is  stored  in  a  3-dimensional  array,  PATCH(48, 
Ml,  M2).  Figure  2  illustrates  the  way  in  which  the  associated  data 
are  stored.  Each  partition  represents  all  the  data  associated  with 
the  corresponding  corner — for  example,  Partition  2  contains  all  the 
data  for  Corner  2. 


Figure  2  -  SMA  Data  Structure 


FUNCTION  LGETV  1 

Entry  Points : 

LGETV,  LSETV 

This  function  retrieves  (or  stores)  surface  patch  information  from 
(or  within)  the  PATCH  data  structure,  given  the  patch  reference  numbers. 


Calling  Sequence: 

N  =  LGETV (PATCH,  Ml,  M2,  IP,  JP,  NPT,  KODE,  VECTOR) 


Parameter  Definitions: 


(Input)  NPT  Number  of  the  corner  to  be  referenced  (Figure  1) . 

If  NPT  =  0,  all  four  corner  points  will  be  referenced. 

KODE  Specific  vector  at  corner  NPT  to  be  referenced. 

0  All  four  vectors 

1  Point 

2  Tangent  along  U-curve 

3  Tangent  along  V-curve 

4  Twist 


(Output)  VECTOR  Array  to  receive  patch  data,  or  that  containing 

information  to  be  stored  into  the  data  structure. 

LGETV  =  0  Function  is  completed 
=  1  Error  has  developed 


Detailed  Description: 

When  LGETV  is  called,  with  NPT  =  0,  the  entire  patch  (48  numbers) 
is  returned  in  VECTOR.  Otherwise,  NPT  indicates  the  specific  corner 
to  be  referenced.  If  all  12  numbers  associated  with  a  corner  are  desired, 
KODE  is  given  as  0.  Otherwise,  KODE  denotes  the  specific  3-dimensional 
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vector  to  be  referenced.  The  vector  obtained  is  positioned  immediately 
following  PATCH(3*(K0DE-1)  +  12*(NPT-1),  IP,  JP)  in  the  Patch  array. 

If  an  attempt  is  made  to  reference  a  surface  patch  outside  of  the 
bounds  of  the  data  structure,  or  if  there  are  any  other  illegal  calling 
parameters,  LGETV  is  returned  with  a  value  of  1.  Otherwise,  its  value 
is  0. 

The  description  of  entry  LSETV  applies  also  for  LGETV,  except  that 
LSETV  updates  data  in  the  structure  whereas  LGETV  retrieves  data  from 
the  structure. 

Future  Modifications; 

In  the  M0D2  version  of  the  SMA  subroutines ,  LGETV  and  LSETV  will  use 
the  Interactive  Data  Manager  for  the  PATCH  data  structure  storage. 

Subroutines  Called:  None. 


FUNCTION  IGETV 

Entry  Points: 

IGETV,  ISETV 

This  function  retrieves  (or  stores)  surface  patch  information  from 
(or  within)  the  PATCH  data  structure,  given  the  curve  numbers  of  an 
intersection. 

Calling  Sequence: 

N  =  IGETV (PATCH,  Ml,  M2,  KNU,  KNV,  NPT,  KODE,  VECTOR) 


Parameter  Definitions: 

(Input)  KNU  U-curve  number  (Figure  1) 

KNV  V-curve  number  (Figure  1) 

NPT  Corner  point  number  to  be  referenced  (Figure  3) 
KODE  Specific  vector  to  be  referenced  (See  LGETV) 

(Output)  VECTOR  Array  to  receive  patch  data,  or  that  containing 
information  to  be  stored  into  the  data  structure. 

IGETV  =  0  Function  is  completed 
=  1  Error  has  developed 


Detailed  Description; 

This  routine  is  similar  to  LGETV;  thus,  only  differences  between 
the  two  will  be  noted. 

(1)  The  U-  and  V-curve  numbers  are  specified  instead  of  the  PATCH 
reference  numbers. 

(2)  NPT  may  not  equal  0 


Subroutines  Called:  LGETV,  LSETV 


1 

- U-CURVE  I 

PATCH  (H,  1-1,  J-l) 

PATCH  (N,  1-1,  J) 

CORNER 

POINT  4 

OANEA 

corner^ 

POINT  2 

POINT  3  ^  u 

^PomT^- '  w 

N— »  V-CURVE  J 

PATCH  (N,  I,  J-l) 

PATCH  (N,  1,  3) 

i 

/ 

■ 

The  corner  point  is  found  in  the  corresponding 

partition  of  each  of  the  four  patches. 

Figure  3  -  A  Shared  Interior  Corner  Point 
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MAIN  COMPUTATIONAL  ROUTINES 


A  mathematical  model  of  a  surface  through  a  given  set  of  points 
can  be  constructed  using  the  two  routines  described  in  this  section. 
There  is  one  condition:  the  set  of  points  must  form  a  quadrilateral 
grid  network.  However,  since  one  or  two  of  the  sides  of  the  patch 
may  be  degenerate,  a  two-  or  three-sided  figure  can  be  obtained.  The 
quantities  determined  in  subroutine  SPLINE  are  endpoint  tangent  rates. 
After  the  endpoint  tangent  rates  have  been  determined,  interior  points 
of  patches  can  be  computed  with  GENSURF. 

The  techniques  used  in  the  construction  of  the  surface  are  described 
in  Appendix  A  which  contains  a  paraphase  of  sections  of  a  paper  by 
I.  M.  Yuille^ 


|  SUBROUTINE  SPLINE  ' 

This  routine  coordinates  all  of  the  operations  needed  to  be  performed 
in  order  to  fit  a  SPLINE  curve  with  data  from  an  array  of  surface  patches. 

Calling  Sequence: 

CALL  SPLINE (PATCH,  Ml,  M2,  CURVPT,  RATES,  ENDS,  SPACE,  ID,  KRV, 

IERR) 


Yuille,  I.  M. ,  "A  System  for  On-Line  Computer  Aided  Design  of 
Ships  -  Prototype  System  and  Future  Possibilities,"  paper  presented 
at  meeting  of  Royal  Institution  of  Naval  Architects,  London  (Apr 
1970) . 
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Parameter  Definitions.* 


CURVPT(J) 

An  array  used  to  store  points  retrieved 
data  structure  (J  <  3*  MAX  (Ml,  M2)) 

from  the 

RATES (J) 

An  array  in  which  parametric  tangents  are  stored 
before  entering  them  into  the  data  structure. 

ENDS (8) 

An  array  used  for  storing  endpoint  tangent  conditions. 

SPACE (K) 

A  scratch  array  used  in  the  solution  of 
equations  in  tangent  calculations.  (K  < 

simultaneous 

4*  MAX (Ml,  M2,  3)) 

ID 

Type  of  Curve, 

ID  =  1  for  U-curve,  ID  =  2  for  V-curve 

KURV 

Curve  Number. 

IERR 

Number  of  retrieval  and  updating  errors  that  occurred 
during  manipulation  of  patch  information.  (Should 
be  0  if  calling  parameters  are  good.) 

Detailed  Description: 

SPLINE  assumes  that  the  patch  array  was  initialized  to  -0  (all  bits  on 
77777777777777777777B)  at  the  programs  start.  The  reason  for  this  initial 
ization  is  that  SPLINE  checks  each  tangent  vector  on  the  curve  KURV  to 
determine  if  any  tangents  were  predefined  by  the  user.  If  tangents  have 
been  defined,  or  if  SPLINE  finds  any  point  discontinuities,  the  curve  is 
broken  into  segments  with  the  defined  tangents  as  endpoint  tangent  con¬ 
ditions.  If  SPLINE  encounters  a  defined  tangent  at  the  beginning  of  an 
interval  but  not  at  the  end  of  the  previous  interval,  the  curve  is  broken 
and  a  free  end  is  assumed  for  that  previous  interval.  The  same  is  done 
for  a  defined  tangent  at  the  end  of  an  interval  and  an  undefined  beginning 
of  the  next  interval  (See  Figure  4) .  If  SPLINE  encounters  point  discon¬ 
tinuities  with  no  defined  tangents  at  the  points,  free  ends  are  assumed. 
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Figure  4a  -  If  there  are  no  defined 
tangents  (except  possibly  the  curve 
endpoints) ,  the  first  derivative 
will  be  continuous. 


'Figure  '4b  -  'ft  a  tangent  is 'defined 
at  the  final  point  of  a  segment,  a 
free  (hinged)  end  is  assumed  for  the 
succeeding  segment  (unless  it,  too, 
is  defined).  In  general,  slope  will 
be  discontinuous. 


Figure  4c  -  If  a  tangent  is  defined 
at  the  initial  end  of  a  segment,  the 
final  end  of  the  preceedlng  segment 
is  assumed  to  have  a  free  (hinged) 
end  (unless  it  is  defined).  Ih 
general,  slope  will  be  discontinuous* 


Figure  4d  -  If  both  tangents  at  a 
single  point  are  defined,  the 
segments  will  have  the  corre¬ 
sponding  slope;  if  the  tangents 
are  scalar  multiples  of  each  other, 
the  slope  will  be  continuous;  if 
the  tangents  have  equal  magnitude 
and  opposite  direction,  the  first 
derivative  will  be  continuous. 


Figure  4  -  Spline  Segments  with  Interior  Tangents  Defined 
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Cross-curve  point  discontinuities  are  handled  by  spline  fitting  the 
curve  twice.  For  example,  if  a  U-curve  was  being  splinefit,  the 
V  =  0  boundary  curve  would  be  fit  first,  then  the  V  =  1  boundary 
curve  would  be  fit.  (Figure  5). 

In  general,  SPLINE  will  accommodate  any  kind  of  surface  or  slope 
discontinuity  as  long  as  the  discontinuity  is  on  a  patch  boundary  or 
at  a  corner  point. 

Subroutines  Called:  CUBIT,  INCONP,  ISETV,  IGETV,  EPTC,  MZCAM 


SUBROUTINE  GENSURF 

This  routine  generates  points  on  a  surface  within  a  particular 
patch,  given  the  surface  tensor  of  the  patch. 

Calling  Sequence: 

CALL  GENSURF  (PATCH,  Ml,  M2,  IP,  JP,  XYZ,  JVSP,  IVSP) 

Parameter  Definitions; 

(Input)  JVSP,  iysp  Number  of  sub-intervals  in  the  V  and 

U  directions,  respectively  +1. 

XYZ (3,  JVSP,  IUSP)  Array  containing  intermediate  boundary 

and  surface  points. 
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Detailed  Description: 

This  routine  computes  each  point  on  the  surface  by  evaluation  of 

t  t 

the  surface  equation:  (X,  Y,  Z)  =  UMBM  V 


where 


U  =  [U3  U2  U  1]  0  <  U  <  1 
V  =  [V3  V2  V  1]  0  £  V  <  1 


■2-2  11 
-3  3-2  1 
0  0  10 
10  0  1 


M  = 


B  = 


00 

01 

00 

V 

01 

V 

10 

11 

10 

V 

11 

V 

00 

u 

01 

u 

00 

uv 

01 

uv 

10 

u 

11 

V 

10 

uv 

11 

w 

Patch  Surface  Tensor 


First,  MBMfc  is  computed.  Then  V  is  varied  from  0  to  1  for  each  value 
of  U  from  0  to  1.  The  increments  in  the  U  and  V  direction  are  1/IUSP 
and  1/JVSP,  respectively.  GENSURF  uses  a  labelled  COMMON/ SCRATCH/  <that 
Is  184  words  in  length  )  for  working  space. 


Subroutines  Called: 
LGETV,  MULT 
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GENERAL  PURPOSE  SUBROUTINES 


|  SUBROUTINE  CUBIT  | 


This  routine  fits  a  parametric  cubic  spline  through  a  series  of 
points  in  an  n-dimensional  space  such  that  specified  end  conditions 
are  satisfied. 


Calling  Sequence: 

CALL  CUBIT (PTS,  RATES,  NDIM,  ENDS,  SPACE,  MAXPTS,  PRMNT) 


Parameter  Definitions: 

(Input)  PTS  (NDIM,  MAXPTS) 

NDIM 

ENDS (J) 

SPACE (K) 

MAXPTS 

PRMNT 

(Output)  RATES  (NDIM,  MAXPTS) 


Array  containing  points  through 
which  a  cubic  spline  will  be 
fitted.  Points  must  be  in 
succession. 

Number  of  dimensions  in  point 
space  (and  first  dimension  in 
the  rates  and  PTS  arrays) . 

Endpoint  conditions  for  the 
spline  fit.  (J  =  2*  NDIM  +2) 

See  note  on  endpoint  condition 
setup. 

Array  used  as  a  working  area  for 
the  solution  of  MAXPTS  simultaneous 
equations .  (K  <  R*  MAXPTS) . 

The  number  of  points  to  be  fitted. 

The  range  of  the  parameter  over  a 
curve  segment  as  it  goes '".from  the 
beginning  point  to  the  end  point. 

Array  containing  endpoint  tangent 
rates. 


Detailed  Description: 

The  SPLINE  approximation  is  a  technique  which  attempts  to  reproduce 
mathematically  a  draftsman's  method  of  drawing  fair  curves. 
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The  differential  equation  for  the  mechanical  spline  curve  is 


where 


Y" 

'2 

(1  +  Y  ; 


3/2 


M(x) 

El 


Y  =  f  (x) 


(1) 


M(x)  Bending  moment  in  the  Spline. 

E  Youngs  modulus. 

I  Cross  sectional  moment. 

In  linerized  beam  theory,  Y*  is  assumed  to  be  small  compared  to  1. 
Therefore,  Equation  1  reduces  to 

Y"  =  f"  (x) 

Clearly,  Y  can  be  expressed  as  a  cubic  polynomial.  Generalizing  this, 

Y  and  X  could  be  expressed  as  cubic  polynomial  functions  of  a  parameter 
U.  This  routine  fits  each  component  of  an  n-dimensional  spatial  set  of 
points  with  a  set  of  parametric  cubic  polynomials  that  is,  for  example, 
xi  =  X±(u)  i  =  1,  2,  .. .  n  -1 

X_t  =  Y±(u)  0  <  U  <  1 

Z±  =  Z±(u) 
for  Cartesian  coordinates. 

A  parametric  spline  function  with  joints  IL  ,  U„,  U„,  ...  U  is  a 

12  3  n 

function  S(U)  constructed  such  that: 

1)  In  each  of  the  Intervals  [U^,  I^],  [t^,  U^] ,  ...  [Un_^-, 

TJn] ,  S(U)  is  a  cubic  polynomial, 

2)  S(U)  has  continuous  first  and  second  derivatives  over  the 
interval  [U, ,  U  ] , 

3)  If  S^Cu)  is  a  spline  approximation  to  x  =  f (u) ,  then 
S  (u)  =  f (u)  at  all  joints  in  the  considered  interval. 
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_  til 

A  cubic  polynomial  (u) ,  (the  i  interval  of  the  spline  function 
til 

for  the  j  dimensional  components  of  the  points),  can  be  found  in  a 
straight-forward  manner  if  ±  <U±) ,  SjPui+P ,  S^pup,  and  ^.(Ui+1)  are 
known,  where  S.pu.)  and  S.pu.+p  are  point  components  and  SY.(U.,p  are 
endpoint  rate  or  tangent  components.  For  "PRMNT"  =1.,  the  following 
results  are  obtained:  (See  Figure  6,  page  22,  for  notation) 


C±1  =  PTS  (j  ,  i) 
c±2  =  RATES  (j ,  i) 

C±3  =  3*  (PTS  (j  ,  i+1)  -  PTS(j,i))  -  2*  RATES  (j  ,i)  -  RATES  (j  ,i+l) 
C±4  =  2.*  (PTS(j,i)  -  PTS(j ,  i+1))  +  RATES  ( j ,  i)  +  RATES (j ,i+l) 


4 

=  Z 
K=1 


C..  H(k  '» 
ik 


i  =  1,  2,  ...  n  -1 
0  <  U  <  1. 


where:  ILpup  «  PTS(j,i) 

^ji(Ui+l)  =  PTS0>i+1> 
S^pU.)  =  RATES  (j,i) 

SjiCb'p  ■  =  RATES 05, i+1) 


(2) 


For  this  program,  the  lengths  of  all  of  the  intervals  [TJ . ,  U.. are 

1  1+1 


equal  to  "PRMNT". 

A  brief  description  of  the  calculation  of  the  entries  of  the  rates 
array  will  now  be  given.  A  necessary  condition  for  S(u)  to  have  con¬ 
tinuous  first  and  second  derivatives  is  that  the  equation  below  must 
be  satisfied. 


19 


s’ (U±  -1)  +  4  s'  (u±)  +  s’(ui+1)  -  |  JS(U1+1)  -  S(U±  _x)] 


where  SL  =  PRMNT  (3) 

2  <  i  <  n  -1 

For  a  problem  with  n  joints,  these  n-2  equations  must  be  supplemented 
with  two  endpoint  tangent  conditions.  If  the  ends  are  free, 

S"(U^)  =  S"(U  )  =  0.  This  condition  is  equivalent  to  the  following 
two  equations : 

4  S'  (u1)  +  s'(u2)  =  |  ls(u2)  -  SCU^l 

4  S '  (U  )  +  S '  (U  .)  -  J  [S(U  )  -  S(U  )]  (4) 

n  n  -i  £  n  n  — i 

If  the  ends  are  fixed,  then  the  values  from  ENDS  array  are  used  as  the 
two  supplementary  conditions. 

The  matrices  associated  with  the  simultaneous  linear  equations  are 
tridiagonal  and  diagonally  dominant.  The  equations  can  be  written  in 
the  following  form: 


°21  S’<V  +C31  S'<V 

cij  s'(Dj  +c2j  s-cu3)  H-  c3J  s'(u.+l) 

j  =  2,  3,  . 


.  n  -  1 


C.  S'  (U  ..)  +  c_  S'  (u  ) 

In  vn_i.  2n  v  n 


C4n  (5) 


The  coefficients  of  Equation  5  correspond  to  those  shown  in  Equations 
3  and  4. 


A  Gaussian  elimination  method  can  be  used  to  solve  the  equations. 
First,  eliminate  the  lower  diagonal  in  the  following  way: 
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Finally,  the  remaining  unknowns  can  be  found  by 

S'  (U±)  =  [C’41  -  C*3i  *  S' (U±+1)]/C,2i  i  -  n . -1 ,  n  -  2....1 

This  equation  yields  all  of  the  S'(U^)*s  or  RATES  known. 

If  PRNMT  $  1,  then  the  RATES  array  is  scaled  by  PRNMT. 

This  procedure  is  carried  out  for  each  dimension  of  PTS  array. 
Subroutines  Called: 

None 


Endpoint  condition  setup  for  CUBIT.  The  number  of  endpoint  conditions 
is  dependent  on  the  number  of  dimensions  in  the  point  space.  For 
example,  for  2-dimensional  space  (x,y),  there  would  be  six  end  con¬ 
ditions.  The  following  table  gives  information  for  3-dimensional 
space  stored  as  (x,y,z)  in  the  PTS  array,  with  the  parameter  U,  in 
the  interval  (a,b). 
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Free  Ends 


Fixed  Ends 


Pa  ENDS  (1)  =  1. 


ENDS (1)  -  2. 


ENDS (3  through  5)  ■  0. 


ENDS(3) 


ENDS (4) 


Pb  ENDS(2)  =  1. 

ENDS (6  through  8)  *  0. 


ENDS (5)  =  g^|Pa 


ENDS (2)  =  2. 


ENDS(6)  -  ~|Pb 
ENDS (7)  =  |^|Pb 


ENDS  (8)  =  ~|Pb 

Note  that  infinite  slopes  can  be  produced  since,  for  example, 

fx/fu  *  fx  and  could  be  zero.  Also  note  that  scalar  multiples 

of  the  endpoint  conditions  will,  in  general,  produce  different  results, 
even  though  the  derivatives  such  as  would  be  the  same. 


Figure  6  -  Notation  used  in  Subroutine  CUBIT 
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FUNCTION  ICONP  j 

Entry  Points : 

ICONP,  IGONT 

ICONP  compares  two  3-dimensional  vectors  to  determine  if  they  are 
identical.  ICONT  determines  whether  the  two  vectors  are  scalar 
multiples  of  each  other. 

Calling  Sequence; 

J  =  ICONP  (Rl,  R2) 

Parameter  Definitions: 

Input:  Rl,  R2  -  Vectors  to  be  compared 

Output:  ICONP  =  0  if  vectors  are  identical 
=  1  if  vectors  are  distinct 
or 

=  0  if  vectors  are  scalar  multiples 
=  1  if  vectors  have  different  directions 

Detailed  Description: 

ICONP  compares  the  like  components  of  the  vectors  to  determine 
whether  they  are  equal. 

ICONT  first  looks  for  any  zero  components  in  R2  that  do  not  have 
corresponding  zero  components  in  Rl.  If  there  are  no  non-corresponding 
zero  components,  a  scale  factor  D  =  R1/R2  is  computed  and  the  rest  of 
the  components  of  Rl  and  R2  are  checked. 

Subroutines  Called: 

None 
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SUBROUTINE  EPTC  | 

This  routine  sets  up  the  endpoint  tangent  conditions  array  for 
CUBIT. 

Calling  Sequence: 

CALL  EPTC  (ENDPT,  VECTOR,  IAORB) 

Parameter  Definitions: 

ENDPT(8)  -  Endpoint  tangent  conditions  array 

VECT0R(3)  -  Endpoint  tangent  vector 
IAORB  -  Beginning  or  end  of  a  segment 

IAORB  =  1  for  beginning 
IAORB  =  2  for  end 

Detailed  Description: 

EPTC  tests  vector  to  determine  if  it  has  been  defined  by  the  user. 
If  vector  has  all  components  of  "-0."  (or  77777777777777777777B) ,  it 
has  not  been  predefined,  and  the  segment  is  assumed  to  have  a  free  end. 
Otherwise,  the  end  is  assumed  to  be  fixed.  EPTC,  then,  sets  up  the 
ENDPT  array  in  accordance  with  the  convention  outlined  in  CUBIT. 

Subroutines  Called: 

MZCAM 
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I  FUNCTION  MZCAM 

This  function  compares  a  three  dimensional  vector  to  "-0."  or 
777 777777777777 777 77B. 

Calling  Sequence; 

J  =  MZCAM  (R) 

Parameter  Definitions: 

Input:  R(3)  -  Vector  to  be  compared 

Output:  MZCAM  =  1  if  R  =  (-0.,  -0.,  -0.). 

=  0  otherwise 

Detailed  Description: 

Since  a  -0.  =  0.  in  the  computer,  this  function  is  needed  to 
distinguish  between  them.  The  routine  compares  R  to  -0.  by  taking 
one  half  a  word  (30  bits)  at  a  time  and  comparing  that  to 
7777777777B. 

Subroutines  Called: 

None 


SUBROUTINE  MULT  ) 

This  routine  performs  matrix  multiplication  of  any  compatible 
arrays . 

Calling  Sequence: 

CALL  MULT  (ARR1,  ARR2,  ARR3,  NR1,  NC1R2 ,  NC2) 
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Parameter  Definitions: 


Input : 


Output : 


ARR1 (NRl ,  NC1R2) 
ARR2(NC1R2,  NC2) 
NRl 
NC1R2 

NC2 

ARR3 (NRl ,  NC2)  = 


-  Left  array  in  multiplication 

-  Right  array  in  multiplication 

-  Number  of  rows  in  ARR1  and  ARR3 

-  Number  of  columns  in  ARR1  and  number 
of  rows  in  ARR2 

-  Number  of  columns  in  ARR2  and  ARR3 
[ARR1]  X  IARR2] 


Subroutines  Called: 
None 


|  SUBROUTINE  BNDRYPT  *| 

This  routine  returns  the  boundary  point  position  vector  at  a 
given  parameter  value. 

Calling  Sequence: 

CALL  BNDRYPT  (PATCH,  Ml,  M2,  IP,  JP,  ISIDE,  UV,  CURVAL) 

Parameter  Definitions: 

Input:  ISIDE  -  The  number  of  the  boundary  on  which  an  inter¬ 

mediate  value  is  desired  (see  Figure  2) 

UV  -  Parameter  value,  0  <  UV  £  1 

Output:  CURVAL  -  Array  to  receive  point  position  vector 

Detailed  Description: 

Subroutine  BNDRYPT  retrieves  point  and  tangent  information  from 
the  data  structure.  The  coefficients  of  a  cubic  polynomial  describing 
the  boundary  are  then  computed.  (Formulas  appear  in  detailed 
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description  of  subroutine  CUBIT) .  The  polynomial  is  then  evaluated 
at  the  parameter  value  "UV". 

Subroutines  Called: 

LGETV 


FUNCTION  FGBLEND  | 

This  function  returns  the  Coons'  Surface  Blending  Function 
Value  at  a  given  parametric  value. 

Calling  Sequence: 

X  =  FGBLEND  (IFG,  101,  C) 

Parameter  Definitions: 

IFG  Type  of  blending  function 

IFG  =  1  for  F-type 
IFG  =  2  for  G-type 

101  Type  of  F  or  G  function 

101  =  0  for  0-type 
101  =  1  for  1-type 

C  Parameter  value  at  which  the  function  will  be  evaluated 

(0  £  C  £  1) 

Detailed  Description: 

FGBLEND  returns  the  value  of  one  of  the  following  functions: 
Fq(C)  =  2*C3  -3*C2  +1. 

F1(C)  =  -2*C3  +  3*C2 

Gq(C)  =  C3  -2*C2  +  C 

G1(C)  =  C3  -C2 

Subroutine  Called:  None 
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APPENDIX  A 


COONS’  SURFACE  PATCH  REPRESENTATION  TECHNIQUE 

The  following  material  was  taken  from  the  paper  by  I.  M.  Yuille 
entitled  "A  System  for  On-Line  Computer  Aided  Design  of  Ships — Prototype 
System  and  Future  Possibilities,"  presented  at  the  meeting  of  the  Royal 
Institution  of  Naval  Architects  in  London  in  April  1970.  The  letter  v 
has  been  substituted  for  the  letter  w  throughout,  to  allow  the  notation 
to  be  consistent  with  the  rest  of  the  report. 

"The  complete  surface  is  represented  by  a  series  of  smaller  surfaces, 
called  patches,  which  join  together  at  their  edges  with  any  desired  order 
of  continuity.  Local  changes  can  therefore  be  made  to  individual  patches 
or  to  the  whole  surface  depending  on  the  extent  and  nature  of  the  con¬ 
tinuity  conditions. 

"Each  patch  is  represented  by  a  three-dimensional  vectof  equation 
in  two  parameters  u  and  v  so  that  a  point  in  the  surface  is  given  by 
x  =  X(u,v) 

y  =  Y(u,v)  0  <  u  <  1 
z  =  Z(u,v)  0  <  v  <  1 

or  in  vector  notation  by  P(x,y,z)  =  P(u,v).  It  is  convenient  to  omit 
the  commas  and  write  P(uv)  or  simply,  (uv) .  Derivatives  will  be  denoted 
by  vector  expressions  such  as 


"The  basis  of  Coons'  work  is  that  each  patch  is  represented  by  an 
equation  of  the  form 

P(uv)  =  f(uv)  +  g(uv)  +  h(uv)  +  ... 

"A  patch  has  four  boundary  curves  denoted  vectorially  by  (uQ)>  (ul) , 
(Ov) ,  (lv)  and  four  corner  points  (00) ,  (01) ,  (10) ,  (11) . 
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"  The  first  kind  of  Coons'  surface  f  (uv)  is  represented  in  matrix 


form  by 

f(uv)  =  -  [-1  Fq(u)  F^u)] 

0  (uO)  (ul) 

"-1 

(Ov)  (00)  (01) 

Vv) 

(lv)  (10)  (11)_ 

Fx(v) 

where  the  functions  Fq(u),  F^(u),  Fq(v),  Fj(v)  are  called  blending' 

functions.  To  ensure  that  the  surface  passes  through  the  -four  boundaries 

and  the  corner  points  the  blending  functions  must  satisfy  the1 conditions 

Fq(0)  =1  F1(0)  =  0 

Fq(1)  -  0  F1(l)  =  1 

The  blending  functions  are  also  chosen  to  satisfy  the  conditions 

F'(°)  =  F' (1)  =  Fj  (0)  =  F'(l)  =  0 
u  u  $ 

so  that  the  derivative  with  respect  to  v  at  v  =  0,  for  example,  is 
f(uO)v  -  (00)vFq(u)  +  (10)vF1(u) 

It  is  seen  that  the  value  of  the  derivative  across  a  boundary  of  the 
patch  depends  only  on  the  tangent  vectors  at  the  two  ends  of  the 
boundary  and  on  the  form  of  the  blending  functions.  Note  that  the 
derivative  is  independent  of  the  shape  of  the- boundary  curves. 

"  The  second  kind  of  Coons'  surface  is  represented  in  matrix  form 
by 


g(uv)  =  -  [-1  Gq(u)  G1(u)] 

0 

g(u0)v 

g(ui)v 

"-1 

g(0v)u 

(00) 

uv 

<01)uv 

G0(v) 

• 

g(lv)u 

(10) 

uv 

<ll)uv 

G;  (v) 

.  1  J 

The  purpose  of  this  surface  is  to  provide  a  correction  to  be  added 
to  a  surface  of  the  first  kind  in  such  a  way  that  the  derivatives 
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across  the  boundaries  may  be  controlled  without  altering  the  shapes  of 
the  boundaries  already  defined.  The  derivative  with  respect  to  v  at 
the  boundary  (uo)  is 

(u0)v  =  f(uO)v  +  g(uO)v 

and  the  elements  such  as  g(uO)^  are  chosen  so  that  for  example 
g(uO)v  =  (u0)v  -  (00) v  Fq(u)  -  (10)v  F1(u) 

In  order  that  the  blending  functions  shall  not  change  the  positions 
of  the  boundaries  they  must  satisfy  the  conditions 
Gq(0)  -  G0(l)  =  G1(0)  =  G1(l)  =  0 

In  order  to  provide  the  required  derivatives  across' the  boundary  they 
must  also  satisfy  the  following  end  conditions 
Gjj(0)  =1  G^(0)  =  0 

GJ(1)  =  0  G'(l)  *  1 

Note  that  the  cross  derivatives  at  the  corners  arise  in  the  definition  of 
this  second  kind  of  surface. 

'In  an  analogous  manner  a  third  function  h(uv)  can  be  defined  which 
gives  control  of  the  second  derivatives  across  patch  boundaries,  and 
so  on. 

"There  is  no  restriction  on  the  shape  of  the  boundaries  of  a  patch 

except  that  they  should  intersect  at  its  corners  and  can  be  represented 

in  parametric  form.  It  is  convenient  to  describe  the  boundaries  in  terms 

of  the  blending  functions.  For  ,-example ,  theboundary  (uO)  is  given  by 

(uO)  =  (00)  F Q(u)  +  (10)  F1(u) 

+  (00)u  GQ(u)  +  (10) u  G1(u) 

"  The  shape  of  the  surface  represented  by  a  patch  is  a  function  of 
the  position  and  the  derivatives  at  the  corner  points  only  and  of  the 
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equations  chosen  for  the  blending  functions.  The  important  feature  of 
this  work  is  that,  if  a  surface  is  represented  by  a  number  of  contiguous 
patches  having  common  values  of  position  and  derivatives  at  the  corner 
points  and  using  similar  blending  function  equations,  there  will  be 
continuity  across  the  patch  boundaries  of  position  and  slope  if  the 
first  and  second  kind  of  blending  functions  are  used,  and  of  second 
derivatives  as  well  if  second  order  continuity  is  maintained  along 
the  patch  boundaries.  In  the  Master  Geometry  programs »  the  unconstrained 
tangent  vectors  are  chosen  by  fitting  spline-  curves  through  the  corner 
points  so  that  second  order  continuity  is  achieved  along  the  boundaries 
at  those  corners  and  across  adjacent  surface  patch  boundaries. 

"in  ship,  aircraft  and  car  design  it  is  required  that  the  surface 
used  should  be  fair.  No  precise  definition  is  available  but  it  is 
generally  agreed  that  there  should  be-  continuity  of  first  and  second 
derivatives  (and  no  unintended  points  of  inflection).  It  has  been 
found  in  practice  that  a  fair  surface  may  be  represented  very  well 
indeed  by  using  blending  functions  of  the  first  and  second  kinds 
only  and  by  taking  the  form  of  these  functions  to  be  cubic  polynomials 
as  follows : 


F  (s)  = 

3  2 

2s  -  3s  + 

o 

F^s)  = 

3  2 

-2s  +  3s 

G  (s)  = 

s3  -  2s2  + 

o 

G1(s)  = 

3  2 

s  -  s 
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"The  equations  may  then  be  combined  and  written  in  the  form 


(uv)  =  [Fq(u)  F1(u)  Gq(u)  G1(u)] 


(00)  (01)  (00), 


(10)  (11)  (10) 


v 


(01), 


(00)  (01)  (00)  (01) 

u  u  uv  t 


(10)  (11)  (10)  (11) 
u  u  uv  uv 


Vv) 

^(v) 

Gq(v) 

Gx(v) 
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This  form  of  Coons'  surface  may  thus  be’ described  by  a  bicubic  para¬ 
metric  equation.  By  including  second  order  continuity  along  the  patch 
boundaries  and  using  the  cross  derivative  terms  fair  surfaces  may  be 
described  for  which  both  the  first  and  second  derivatives  are  continuous 
throughout.  In  general  this  is  necessary  for  the  satisfactory  represen¬ 
tation  of  doubly  curved  surfaces.  Note  that  the  coefficients  in  the 
square  matrix  include  derivatives  with  respect  to  the  parameters  u  and 
v  and  do  not  represent  physical  tangents  at  points  on  the  surface.  The 
latter  may  easily  be  obtained,  however,  by  means  of  equations  of  the  type 

dy  _  dj_  du  _9y.  jjv 
dx  9u  9x  8v  8x 

Surfaces  described  in  this  way  may  have  virtually  any  shape  and  the  method 
is  eminently  suited  to  engineering  design.  Surface  description  by  para¬ 
metric  equations  has  several  advantages  compared  with  non-parametric 
description.  In  particular  the  equations  are  axis-independent,  so  that 
transformations  of  axes  or  scale  are  easily  computed,  and  avoid  the 
possibility  of  infinite  slopes  such  as  may  occur  if  non-parametric 
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equations  are  used.  The  analytical  and  arithmetical  techniques  used, 
for  example  to  calculate  intersections,  are  rather  more  complex  than 
with  other  methods.  The  technique  is  intended  for  use  on  a  computer, 
however,  and  once  it  has  been  implemented  the  designer  is  relieved 
of  the  necessity  to  follow  the  analytical  geometry  in  detail.  .  . " 
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APPENDIX  B 


TEST  RUNS  OF  SMA  SUBROUTINES 

The  performance  of  the  SMA  routines  was  tested  and  verified  by  an 
application  program.  This  program  consists  mainly  of  calls  to  SPLINE, 
GENSURF,  and  THREED  (a  routine  to  do  3-dimensional  plotting  on  the 
CalComp  765  plotter).  Figures  7,  8,  and  9  are  output  of  this  test 


program. 
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Figure  9  -  An  Arbitrary  25-Patch  Surface  Produced  by  SMA  (Example  B) 
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